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We  present  two  experimental  schemes  that  can  be  used  to  implement  the  Factorized  Quantum 
Lattice-Gas  Algorithm  for  the  ID  Diffusion  Equation  with  Persistent-Current  Qubits.  One  scheme 
involves  biasing  the  PC  Qubit  at  multiple  flux  bias  points  throughout  the  course  of  the  algorithm. 

An  implementation  analogous  to  that  done  in  Nuclear  Magnetic  Resonance  Quantum  Computing 
is  also  developed.  Errors  due  to  a  few  key  approximations  utilized  are  discussed  and  differences 
between  the  PC  Qubit  and  NMR  systems  are  highlighted. 
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I.  INTRODUCTION 

Most  algorithms  designed  for  quantum  computers  will  not  best  their  classical  counterparts  until  they  are  imple¬ 
mented  with  thousands  of  qubits.  For  example,  the  factoring  of  binary  numbers  with  a  quantum  computer  is  estimated 
to  be  faster  than  a  classical  computer  only  when  the  length  of  the  number  is  greater  than  about  500  digital.  Account¬ 
ing  for  error  correction  circuitry!  would  bring  the  size  of  the  needed  quantum  computer  to  be  in  the  thousands  of 
qubits.  In  contrast,  the  Factorized  Quantum  Lattice-Gas  Algorithm  (FQLGA)!  for  fluid  dynamics  simulation,  even 
when  run  on  a  quantum  computer  significantly  smaller  than  the  one  just  discussed,  has  significant  advantages  over 
its  classical  counterparts. 

The  FQLGA  is  the  quantum  version  of  classical  lattice-gases  (CLG)!.  CLG  are  an  extension  of  classical  cellular 
automata  with  the  goal  of  simulating  fluid  dynamics  without  reference  to  specific  microscopic  interactions.  The 
binary  nature  of  the  CLG  lattice  variables  is  replaced  for  the  FQLGA  by  the  Hilbert  space  of  a  two-level  quantum 
system.  The  results  of  this  replacement  are  similar  to  that  of  the  lattice-Boltzmann  model,  but  with  a  few  significant 
differences!.  The  first  is  the  exponential  decrease  in  required  memory.  The  second  is  the  ability  to  simulate  arbitrarily 
small  viscosities. 

As  of  today  there  is  a  plethora  of  qubits  to  choose  from  when  designing  a  quantum  computer,  and  a  promising 
class  is  superconducting  qubits  based  on  Josepshon  junction  circuits6'7'8'9,10.  One  major  advantage  of  any  of  these 
superconducting  systems  is  the  ability  to  precisely  engineer  the  quantum  Hamiltonian,  which  extends  from  single 
qubit  design  to  multi-qubit  coupling  arrangements  to  measurement  engineering.  The  quantum  computer  considered 
here  will  be  built  using  the  Persistent-Current  Qubit  (PC  Qubit)!. 

The  goal  of  this  paper  is  to  show  how  one  can  implement  a  one  dimensional  version  of  the  FQLGA  with  the  PC  Qubit. 
To  this  end  we  will  begin  by  reviewing  the  algorithm,  specifically  the  one  that  simulates  the  diffusion  equation,  without 
a  loss  of  generality  in  understanding  the  essence  of  the  algorithm  or  its  general  requirements.  We  will  then  review  the 
PC  qubit  and  show  explicitly  how  to  implement  the  algorithm  with  this  system.  Some  important  differences  between 
the  PC  qubit  and  the  two-state  system  studied  in  Nuclear  Magnetic  Resonance  Quantum  Computation  (NMRQC)!! 
will  be  shown  to  allow  for  some  interesting  new  techniques  in  implementing  quantum  logic.  We  will  also  show  how 
to  implement  the  algorithm  with  the  PC  qubit  in  a  very  analogous  way  to  NMRQC  schemes!!,  with  a  few  significant 
differences. 


II.  FQLGA  FOR  THE  ID  DIFFUSION  EQUATION 

The  first  thing  one  must  do  in  the  FQLGA  is  to  define  a  lattice.  Each  lattice  point  n  will  represent  a  unique  position 
in  the  simulated  fluid.  The  simulation  will  contain  a  finite  number  of  lattice  points,  hence  space  is  discretized  in  the 
simulation. 
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Next  one  must  encode  the  mass  density  p  of  the  fluid  at  each  lattice  site.  In  the  FQLGA  this  is  done  by  building  at 
each  lattice  site  a  set  {«}  of  coupled  qubits.  Each  qubit  represents  the  motion  of  particles  on  the  microscopic  level  in 
one  of  a  finite  set  of  directions.  For  the  diffusion  equation  in  one  dimension,  at  any  point  in  your  fluid,  there  are  only 
two  possible  directions  for  each  particle  to  be  moving,  to  the  left  and  to  the  right.  Hence,  only  two  qubits  are  needed 
to  specify  the  mass  density  pn  at  each  lattice  site.  This  intuitive  reasoning  does  not  extrapolate  to  higher  dimensional 
simulations  because  even  in  two  dimensions  there  would  be  an  infinite  number  of  directions  particles  could  travel  in. 
In  higher  dimensions  one  must  adhere  to  much  more  mathematical  conditions  to  decide  on  the  small  set  of  directions 
one  must  include  for  a  faithful  simulation^.  The  probability  P  of  a  particle  to  be  participating  in  the  motion  assigned 
to  each  qubit  will  be  encoded  in  the  probability  amplitude  of  the  qubit  being  in  its  excited  state  |1).  The  state  of  a 
qubit  is  thus  set  to 


iT!n)  =  yr^wrio)  +  v/^rii>  (i) 

where  i  is  the  qubit  index,  n  is  the  lattice  site  index,  and  |0)  is  the  ground  state  of  the  qubit.  For  the  one  dimensional 
problem  considered  here,  *  =  {1,2}  and  n  =  { 1,  N}  where  N  is  the  number  of  lattice  sites  used  in  the  simulation.  One 
can  easily  conceive  of  fluids  of  multiple  phases  with  multiple  types  of  interactions  even  in  one  dimension,  in  which 
the  size  of  {*}  would  be  much  larger,  but  this  will  not  be  considered  here.  The  mass  density  p  is  then  calculated 
by  summing  the  occupation  probabilities  for  all  qubits  at  a  node.  At  time  t=0  in  a  ID  simulation  the  occupation 
probabilities  P\n  and  P2n  are  set  to  pn / 2,  which  is  the  condition  for  local  equilibrium  in  the  fluidic. 

Now  that  the  fluid  is  initialized,  one  must  account  for  the  interaction  of  particles  in  the  fluid.  These  collisions 
are  encoded  by  the  application  of  a  unitary  transformation  to  the  coupled  systems  at  each  lattice  site.  For  the  ID 
diffusion  equation  this  unitary  transformation  is 
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The  basis  for  computation  is  the  set  of  four  product  states:  |0)|0),  representing  no  particles  at  the  site,  |0)|1), 
representing  the  existence  of  only  a  particle  moving  to  the  right  at  the  site,  1 1}  1 0) ,  representing  the  existence  of  only 
a  particle  moving  to  the  left  at  the  site,  and  1 1}  1 1) ,  representing  particles  moving  in  both  directions  at  the  site.  To 
conserve  particle  number  there  can  be  coupling  only  between  the  middle  two  states.  The  identity  transformation  on 
the  first  and  last  states  corresponds  to  no  collisions  and  a  perfectly  elastic  collision  respectively.  Transformation  of 
the  middle  two  states  was  something  that  never  existed  in  the  classical  algorithm  because  there  was  no  superposition 
of  these  two  states. 

After  collision  the  states  of  the  qubits  at  each  lattice  site  are  in  general  entangled,  and  we  denote  that  state  as 
|Yn).  The  state  of  each  qubit  is  then  measured,  denoted  by  |Xi"),  and  the  process  described  thus  far  is  repeated  many 
times  to  achieve  an  ensemble  average.  Upon  completion  of  these  measurements  one  will  have  found  the  post-collision 
outgoing  occupation  probabilities,  denoted  by  Ptn  once  again.  Note  that  the  occupation  probabilities  now  represent 
something  very  different  than  before  the  collision.  The  particles  have  now  interacted  and  are  ready  to  move  to  the 
next  lattice  site. 

One  must  now  “stream”  the  occupation  probabilities  to  their  new  lattice  sites.  This  is  done  in  a  classical  computer 
by  storing  the  occupation  probabilities  at  each  lattice  site  that  are  coming  from  adjacent  lattice  sites  due  to  colli¬ 
sions.  More  precisely,  P\n  becomes  Pi"+1  and  P2n  becomes  P2n  l-  Periodic  boundary  conditions  are  assumed  when 
streaming  at  the  edges  of  the  fluid. 

To  find  the  mass  density  pn  at  t=l  one  simply  adds  the  occupation  probability  for  both  qubits  at  site  n  once 
streaming  has  been  done.  One  time  step  of  the  algorithm  has  now  been  completed.  To  simulate  the  next  time 
step  simply  start  the  above  procedure  all  over  again  except  now  setting  the  initial  states  with  the  new  occupation 
probabilities  just  found. 

The  algorithm  can  be  summarized  by  four  major  steps,  which  are  illustrated  in  figure  ^  The  first  step  encodes 
the  initial  state  of  the  fluid  by  quantum  mechanically  setting  the  state  |H/jn)  of  each  qubit  at  each  lattice  site.  The 
second  step  transforms  the  two-qubit  product  state  at  each  lattice  site  to  in  general  an  entangled  state,  whose  state  is 
denoted  by  |T”).  Third  one  makes  a  projective  measurement  of  the  post-collision  states  | Xin)>  and  one  must  repeat 
the  first  three  steps  to  find  the  outgoing  occupation  probabilities  P™ .  In  the  fourth  and  final  step  one  streams  the 
mass  density  with  the  appropriate  post-collision  occupation  probabilities,  from  the  left  with  particles  representing 
positive  momentum,  and  from  the  right  with  particles  representing  negative  momentum,  and  the  mass  density  is 
calculated.  Subsequent  time  steps  are  identical  except  for  a  change  in  the  initial  mass  density  profile,  i.e. ,  initial  qubit 
states  in  the  first  step. 
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FIG.  1:  General  summary  of  the  four  major  steps  that  comprise  one  time  step  of  the  ID  FQLGA  fluid  dynamics  simulation. 
The  sequence  of  initialization  of  mass  density,  collision  of  particles,  and  measurement  of  post-collision  states  is  repeated  many 
times  to  make  an  ensemble  measurement.  Propagation  between  collisions  is  accomplished  by  storing  the  adjacent  occupation 
probabilities  for  a  given  site  in  a  classical  computer,  where  the  mass  density  is  then  calculated  for  this  time  step.  Subsequent 
time  steps  utilize  these  “streamed”  occupations  when  initializing  again  for  the  next  set  of  collisions  and  “streaming” . 


III.  PERSISTENT-CURRENT  QUBIT 


The  fundamental  unit  of  quantum  logic  we  will  use  to  implement  the  algorithm  is  the  PC  Qubiti4.  It  consists  of  a 
superconducting  loop  that  is  interrupted  by  three  Josephson  junctions,  pictured  as  x’s  in  figure  |5{a).  The  magnetic 
flux  $  is  the  only  control  field  for  our  qubit,  and  as  shown  in  the  figure,  is  usually  denoted  by  /  =  where 

=  h/ 2e  is  a  single  flux  quantum,  h  is  Planck’s  constant,  and  e  is  the  magnitude  of  the  charge  of  an  electron. 
Physically,  a  Josephson  junction  is  a  small  layer  of  insulator  sandwiched  between  superconductors,  so  our  system  is 
a  superconducting  loop  interrupted  by  three  layers  of  insulator  about  lnm  thick.  For  single  qubit  manipulation  the 
magnetic  flux  through  the  loop  will  be  modified.  The  flux  seen  by  a  DC  SQUID  magnetometer,  a  combination  of 
applied  flux  and  qubit-induced  flux,  will  serve  as  our  measurement  variable. 

The  Hamiltonian  of  the  qubit  is  derived  by  considering  a  circuit  element  model  of  our  system,  which  consists  of 
three  Josepshon  junctions,  where  two  junctions  have  the  same  cross-sectional  area,  and  the  third  is  smaller  by  a  factor 
of  a.  The  constituent  relations  for  an  ideal  Josephson  junction  are 


I  =  4  sin(^) 
v  = 

2tt  dt 


(3a) 

(3b) 


where  /  is  the  current  through  the  junction,  V  is  the  voltage  across  the  junction,  Ic  is  the  maximum  current  the 
junction  can  hold  without  a  voltage  appearing  across  it,  ip=6i-02,  and  6*12  is  the  phase  of  the  plane  wave  macroscopic 
wavefunction  that  characterizes  the  superconductor  condensate  on  the  +,-  side  of  the  junction  respectively.  Note  that 
Ic  is  a  function  linear  in  the  cross-sectional  area  of  the  junction,  and  hence  the  third  junction  has  a  lower  Ic  by  a 
factor  of  a. 

The  energy  associated  with  an  ideal  Josephson  junction  is  found  by  integrating  the  power  from  time  t  =  0  to  some 
final  time  t0,  which  is  equivalent  to  an  integral  from  zero  phase  to  some  phase  ip.  The  energy  it  takes  to  set  the  phase 
of  a  Josephson  junction  to  (p  is 


P  fto  n  •  ,w$o  dp'  $0/c 

E=  /  (Ic  sm  <£>)(— — —  )dt  = 


'  2ir  dt 


2tt 


sin  ip'd(j>'  =  Ej(  1  —  cos<p) 


(4) 


where  Ej  =  $>0Ic/2n. 
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FIG.  2:  (a)  Schematic  drawing  of  the  PC  Qubit.  The  x’s  represent  Josephson  junctions,  with  all  connecting  leads  made  of 
the  same  superconductor  that  is  part  of  the  junctions.  The  sign  conventions  chosen  when  summing  phases  are  shown,  and  the 
magnetic  flux  penetrating  the  loop  (in  units  of  $0)  is  labeled  by  /.  (b)  The  potential  energy  of  the  full  Hamiltonian  for  the 
PC  Qubit  is  plotted  when  the  system  is  biased  at  /  =  0.4954>o.  The  phase  particle  sees  an  infinite  2D  lattice  with  unit  cells 
resembling  a  double  well  potential. 


By  including  the  charging  energies  due  to  the  capacitance  of  the  junctions,  the  Hamiltonian  of  our  circuit  isi4 

p2  p2 

H  =  2 \M~  +  2 M~  +  Ej^2  +  a  ~  2cos(^p)  C0S(Pm)  -acos(2nf  +  2 ipm)\  (5) 

where  (pp  =  ^  +  ip2  ,  <pm  =  -  P2,  Pp  =  Mpdipp/dt ,  Pm  =  Mmdipm/dt ,  Mp  =  ($0/2tt)22C,  and  Mm  = 

(<l>0/27r)22C'(l  +  2a).  The  number  of  degrees  of  freedom  in  the  problem  was  reduced  by  the  fluxoid  quantization 
condition^ 


<pi  +  ^3  -  y>2  =  2,-kii  + 


27r$ 


(6) 


which  forces  the  sum  of  the  gauge  invariant  phases  to  be  proportional  to  the  amount  of  flux  quanta  modulo  an  integer 
multiple  of  2-7T. 

We  have  chosen  to  associate  the  capacitive  energy,  the  first  two  terms  in  @,  with  the  kinetic  energy,  and  the  ideal 
Josephson  energy,  the  last  four  terms  in  ®>  with  the  potential  energy.  The  potential  energy  is  that  of  an  infinite 
lattice  of  double  wells,  as  seen  in  figure  0b).  The  arrow  in  the  plot  shows  the  direction  one  would  take  to  traverse 
from  one  side  of  a  double  well  to  another.  The  barrier  between  the  left  and  right  sides  of  a  single  double  well  is  much 
lower  than  any  barrier  between  different  double  wells. 

Though  quantum  mechanics  plays  a  foremost  role  in  deriving  the  constitutive  relations  for  the  superconducting 
circuit  elements,  the  Hamiltonian  for  the  circuit  itself  so  far  has  been  classical.  The  quantum  version  of  the  circuit 
can  be  understood  by  imagining  a  phase  “particle”  in  the  potential  shown  in  figure  0b).  The  behavior  of  this 
“particle”  is  analogous  to  a  particle  with  an  anisotropic  mass  moving  in  a  2D  periodic  potential,  and  so  there  exist 
energy  bands  in  a  fc-space,  which  is  here  related  to  the  charge  stored  capacitively  by  the  Josephson  junctions.  By 
properly  choosing  Ej/Ec,  where  Ec  =  e2/2 C,  one  can  remove  any  k  (and  hence  charge)  dependence  in  the  energy 
of  the  system,  and  hence  can  reduce  the  problem  to  that  of  an  effective  double  well.  What  we  have  done  is  choose 
parameters  such  that  tunneling  between  adjacent  double  wells  can  be  neglected  relative  to  the  tunneling  within  a 
double  well  in  the  tight-binding  solution,  making  the  solution  effectively  that  of  a  single  double  well. 

By  considering  only  the  lowest  two  levels  of  the  double  well,  the  equivalent  Hamiltonian  is 

H  =  $0Ip(f-  ^)dz  -  rdx  (7) 

where  ±Jp  are  the  eigenvalues  of  circulating  current  for  the  two  az  eigenstates  and  r  is  the  tunneling  element  from 
one  side  of  the  double  well  to  the  other.  The  energies  of  the  two  eigenstates  along  with  a  sketch  of  the  double  well 
as  a  function  of  applied  flux  are  shown  in  figure  0  One  significant  difference  between  this  qubit  and  the  one  used  in 
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FIG.  3:  The  energy  levels  of  the  PC  Qubit  are  shown  as  a  function  of  /.  The  eigenstates  of  the  system  change  with  /  and 
are  labeled  on  the  plot.  The  change  in  the  potential  of  the  phase  particle  is  also  depicted  at  the  top  of  the  plot.  The  energy 
difference  between  the  states  at  /  =  1/2  is  seen  to  be  twice  the  intra-well  tunnelling. 


NMRQC  is  the  presence  of  the  ax  term.  The  implication  of  such  a  term  is  that  the  energy  of  the  eigenstates  as  well  as 
the  eigenstates  themselves  change  as  the  bias  field  is  modified.  In  figure|3we  see  that  at  the  classical  degeneracy  point 
/  =  1/2  the  qubit’s  eigenstates  are  ax  eigenstates,  while  far  from  /  =  1/2,  but  still  far  from  /  =  1,  the  eigenstates  are 
those  of  az.  The  same  thing  happens  for  /  <  1/2,  but  now  the  eigenstates  have  switched  energies,  i.e.,  the  ground 
state  here  is  the  first  excited  state  of  az  and  vice  versa. 

The  PC  Qubit  has  some  advantages  over  other  superconducting  qubits.  Charge  fluctuations,  a  consequence  of 
trapped  substrate  charge,  are  deemed  inconsequential  through  the  choice  of  parameters  used  when  designing  the  PC 
Qubit  circuit.  Also,  flux  noise  has  been  reduced  in  this  system  over  other  flux  qubits  since  this  system  has  a  smaller 
loop. 

A  typical  conceptual  misconception  can  be  addressed  at  this  point.  The  two  different  states  used  in  computation  are 
not  related  to  single  Cooper  pair  behavior.  Rather,  they  are  macroscopically  distinct  states  described  by  the  circulating 
current  due  to  millions  of  Cooper  pairs,  characterized  by  different  average  induced  fluxes  when  in  a  magnetic  field. 

As  seen  in  section^  the  qubits  will  need  to  be  coupled.  For  the  PC  Qubit,  just  as  microwaves  can  only  be  coupled 
in  through  az ,  coupling  between  qubits  can  only  be  of  the  form  azaz.  Other  coupling  terms  can  be  introduced  by 
design,  even  for  these  planar  devices^. 

IV.  IMPLEMENTATION  WITH  THE  PC  QUBIT 

We  now  show  how  one  can  use  the  PC  Qubit  to  simulate  the  ID  diffusion  equation.  In  section  HV  Al  we  elaborate 
on  a  scheme  based  upon  changing  the  flux  bias  points  of  the  qubits  during  the  algorithm,  which  will  lead  to  a  very 
general  initialization  scheme,  but  a  less  general  collision.  In  sectionllV  Blwe  discuss  a  more  general  collision,  analogous 
to  that  done  in  NMRQC,  and  how  to  initialize  the  qubits  before  this  general  collision. 


A.  The  Multiple  Bias  Point  Implementation 

The  first  of  the  four  steps  of  the  algorithm  is  initializing  each  qubit  at  each  node.  As  discussed  in  section  [H]  each 
qubit  must  be  initialized  into  a  state  of  real  and  positive  phase  in  its  own  Hilbert  space.  This  set  of  states  consists 
of  all  those  lying  on  the  real  phase  geodesic  between  the  ground  and  first  excited  states  on  the  Bloch  sphere.  The 
ground  state  of  the  PC  Qubit  as  a  function  of  applied  flux  coincidentally  also  occupies  exactly  this  geodesic  on  the 
Bloch  sphere,  as  discussed  in  section  ITTTI  Initialization  can  thus  be  accomplished  while  staying  in  the  ground  state  by 
adiabatically  changing  the  applied  magnetic  flux,  as  depicted  in  figure  |3J 
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Flux  Seen  By  Qubit 


FIG.  4:  The  overlap  between  the  first  (second)  excited  state  |1)  (|2))  of  the  PC  Qubit  coupled  system  and  the  |01)  ( 1 10) ) 
computational  state  are  plotted  when  qubit  1  is  biased  at  f  —  0.508.  The  coupling  between  |1)  and  |2)  in  the  presence  of  an 
AC  magnetic  field  is  also  plotted.  Qubit  coupling  equal  to  r  (same  for  both  qubits)  was  assumed  in  the  calculation. 


The  flux  used  to  set  the  state  of  one  qubit  will  be  affected  by  the  state  of  the  other  qubit  and  its  bias  current. 
This  permanent  inductive  coupling  can  be  accounted  for  by  slightly  adjusting  the  applied  flux  to  compensate  for  the 
flux  introduced  by  the  other  qubit  and  its  bias  line.  All  errors  due  to  approximations  made  when  initializing  by 
rotating  dynamically  on  the  Bloch  sphere,  including  errors  due  to  decoupling,  have  been  avoided^.  We  emphasize 
that  the  initialization  portion  of  the  algorithm  is  identical  for  any  simulation,  whether  it  be  for  a  different  equation, 
a  multi-phase  simulation,  or  in  a  different  number  of  dimensions. 

The  second  step  of  the  algorithm  is  the  collision.  Here  we  study  a  very  specific  unitary  transformation,  the  y/swap 
described  in  section  [H]  This  matrix  simply  “half-way”  swaps  the  middle  two  (first  and  second  excited)  computational 
states  of  the  coupled  system.  In  NMRQC,  the  coupled  eigenstates  are  exactly  those  computational  states,  but  there 
are  no  direct  matrix  elements  connecting  these  states^.  When  the  PC  Qubits  are  coupled,  the  first  and  second  excited 
states  of  the  four- level  system,  denoted  as  |1)  and  |2)  respectively,  are  in  general  not  the  same  as  the  computational 
basis  states  the  swap  intends  to  affect.  However,  the  dc  bias  fields  of  each  qubit  can  be  tuned  to  make  these  two 
sets  of  eigenstates  coincide.  Once  this  is  done,  one  can  then  implement  the  y/swap  by  simply  oscillating  the  magnetic 
field  bias  at  the  frequency  corresponding  to  the  energy  difference  between  the  middle  two  eigenstates.  This  is  just 
a  Rabi  oscillation  between  the  middle  two  eigenstates,  and  since  one  wants  to  only  “half-way”  swap  the  states,  the 
radiation  should  only  be  left  on  for  a  quarter  of  a  Rabi  period. 

Besides  finding  the  appropriate  bias  points  such  that  the  middle  two  eigenstates  of  the  coupled  system  are  very 
similar  to  the  middle  two  computational  states,  one  must  also  verify  that  the  coupling  between  these  states  in  the 
presence  of  an  oscillating  magnetic  field  is  non-zero.  The  results  of  these  calculations  are  shown  in  figure 0]  The  bias 
point  of  qubit  1,  /i,  must  be  chosen  to  be  far  from  1/2,  but  not  too  far.  In  these  calculations  we  take  /i  =  0.508.  In 
the  figure  we  see  that  when  qubit  2  is  biased  at  around  f-2  =  0.51,  the  first  two  system  excited  states  are  very  similar 
to  the  middle  two  computational  states,  with  overlap  elements  of  about  0.97.  At  this  same  bias  point  one  sees  a  Rabi 
matrix  element  of  about  0.02,  which  is  more  than  sufficient  for  our  purposes. 

This  approximate  swap  has  been  incorporated  into  simulation  of  the  FQLGA  for  the  ID  diffusion  equation  and 
the  results  are  pictured  in  figure  O  Snapshots  of  three  different  times  have  been  shown,  for  both  an  ideal  simulation 
and  one  including  the  error  introduced  due  to  the  approximate  collision.  At  time  t  =  0  one  can  see  that  we  have 
initialized  our  fluid  to  a  gaussian  profile.  Later  time  steps  of  the  ideal  implementation  show  the  expected  spreading 
due  to  diffusion,  while  conserving  the  total  number  of  particles.  Increase  in  the  diffusion  constant  of  the  approximate 
collision  when  compared  to  the  ideal  simulation  results  from  the  enhanced  population  in  the  |00),  1 01) ,  and  1 10)  states 
relative  to  the  1 1 1)  state  due  to  extra  matrix  elements  in  the  approximate  swap  that  couple  the  four  states. 

Even  with  an  ideal  swap  operator,  an  interesting  timing  issue  arises  upon  non-adiabatically  switching  the  bias 
fields  from  the  initialization  settings  to  the  proper  settings  to  do  a  Rabi  oscillation  between  the  two  middle  product 
states.  We  first  illustrate  this  timing  issue  and  then  show  how  it  can  be  made  negligible  by  making  a  larger  ensemble 
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FIG.  5:  The  results  of  the  FQLGA  are  simulated  having  accounted  for  the  approximate  nature  of  the  collision  proposed  in 
section  liV  Al-i-l.  Qubit  coupling  equal  to  r  (same  for  both  qubits)  was  assumed  in  the  approximate  swap  simulation.  The  ideal 
results  of  the  FQLGA  are  also  shown(o). 


measurement. 

Once  the  applied  fluxes  are  changed  to  those  appropriate  to  perform  the  approximate  swap,  the  initialized  states 
will  most  likely  not  be  eigenstates  anymore,  and  hence  will  begin  to  precess  due  to  a  time- independent  perturbation. 
Assuming  things  can  not  be  accurately  controlled  at  these  timescales,  one  will  have  now  introduced  a  random  phase 
difference  between  the  two  qubits  due  to  this  Larmor  precession.  This  effect  is  pictured  in  figure  0  The  states  before 
the  bias  fields  are  switched  lie  along  the  same  geodesic.  Upon  changing  the  magnetic  flux  seen  by  each  qubit,  the 
qubits  begin  to  precess,  out  of  phase. 

The  effect  of  this  phase  difference  <5  on  the  algorithm  will  be  to  alter  the  fraction  of  particles  at  each  lattice  site, 
post-collision,  that  are  “moving”  to  the  right  and  to  the  left.  The  results  of  measuring  the  post-collision  occupation 
probabilities  having  accounted  for  a  constant  phase  difference  is  summarized  by 


Pi  =  Pi,  5=0  +  7  sin(<5)  (8a) 

Pi  =  Pi,  s=o  -  7  sin(c5).  (8b) 

The  effect  of  this  error  on  the  simulation  is  effectively  averaged  away  when  an  ensemble  is  measured,  since  S  is  randomly 
different  for  each  member  of  the  ensemble.  These  results  are  shown  in  figure  0  One  can  see  small  random  deviations 
from  the  ideal  simulation  that  can  be  made  infinitesimally  small  by  measuring  a  larger  ensemble  (an  ensemble  average 
of  1000  repeated  measurements  was  simulated  here). 

In  summary,  an  initialization  scheme  has  been  developed  that  is  not  available  to  qubits  with  only  one  term  in 
their  Hamiltonian.  This  initialization  scheme  is  limited  only  by  the  precision  of  the  current  source  used  to  create 
the  magnetic  field  that  biases  the  qubit.  The  scalability  of  this  scheme  relative  to  those  used  in  NMRQC  is  an 
interesting  question,  but  is  not  resolved  here.  The  collision  implementation  is  also  unique  to  qubits  with  multiple 
term  Hamiltonians,  but  the  unitary  transform  implemented  is  unique  to  the  diffusion  equation,  and  fortuitously 
simple.  A  collision  scheme  that  could  be  generalized  to  any  unitary  transformation  would  be  much  more  useful. 


B.  Generalized  NMRQC-like  Implementation 

Generalization  of  the  above  implementation  to  any  fluid  dynamics,  i.e.,  any  unitary  transformation,  can  be  done  in 
an  analogous  way  to  NMRQC  schemes.  Generalization  of  the  collision  transformation  consists  of  using  a  universal  set 
of  quantum  computation  gates,  and  decomposing  all  transformations  into  a  sequence  of  these^.  In  NMRQC  collision 
is  performed  by  a  sequence  of  single  qubit  unitary  transformations  and  coupled  free  evolution.  In  this  section  we  will 
begin  by  discussing  the  single  qubit  rotations  needed  for  a  general  decomposition,  and  briefly  mention  the  role  they 


FIG.  6:  The  unfilled  circle  and  triangle  represent  two  typical  initialized  states  at  one  lattice  point,  both  on  the  same  north 
pole  to  south  pole  geodesic,  before  their  flux  bias  is  changed  to  perform  the  collision.  The  filled  circle  and  triangle  represent 
the  same  states  after  imprecise  bias  changing  has  occurred.  Imprecisely  timed  Larmor  precession  introduces  a  random  phase 
difference  5  between  the  two  states. 
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FIG.  7:  The  results  of  the  FQLGA  are  simulated  having  accounted  for  a  random  phase  difference  introduced  before  the  collision 
for  each  member  of  the  measurement  ensemble(+).  The  ideal  results  of  the  FQLGA  are  also  shown(o). 


could  play  in  initialization.  We  will  then  explore  the  free  evolution  of  a  coupled  PC  Qubit  system,  and  then  show 
how  to  combine  the  single  and  coupled  pulses  to  implement  the  collision  of  the  ID  FQLGA  for  the  diffusion  equation. 

Single  qubit  transformations  can  most  easily  be  achieved  in  a  rotating  frame,  since  here  the  frequency  of  precession 
can  be  much  lower  than  the  Larmor  time  scale.  For  this  implementation  we  will  only  study  the  case  where  our  qubit 
is  biased  at  /  =  1/2.  This  discussion  is  easily  generalized  to  any  bias  point,  but  the  mathematical  notation  can  get 
quite  cumbersome.  The  Hamiltonian  of  the  PC  Qubit  in  an  applied  AC  field  is 

H  =  w0Ix  +  go  cos(w0t  +  <t>)Iz  (9) 

where  w0  is  the  frequency  of  the  applied  field,  g0  is  proportional  to  the  amplitude  of  the  applied  field,  (f>  is  the  phase 
of  the  applied  field,  and  I*  =  hat/2. 
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FIG.  8:  The  results  of  the  FQLGA  are  simulated  for  NMR-like  single  qubit  pulse  initialization,  where  errors  arise  from 
initializing  a  coupled  ground  state  that  is  not  a  product  state(+).  The  ideal  results  of  the  FQLGA  are  also  shown(o). 


In  the  frame  rotating  about  Ix  this  Hamiltonian  becomes 


H  =  -g0[cos((j))Iz  +sin  {4>)Iy\. 


(10) 


The  quantum  state  will  now  precess  in  this  frame  about  the  axis  defined  by  <f> ,  with  the  angle  through  which  the  state 
has  precessed  given  by  9  =  g0t/2.  It  will  be  convenient  to  only  consider  the  set  of  two  rotations  defined  by  (j)  =  0  and 
7t/2,  which  are  rotations  about  I~  and  Iy  respectively.  These  rotations,  denoted  as  Rz(6)  and  Ry{8)  respectively,  can 
be  used  in  conjunction  to  bring  the  qubit  state  to  any  point  on  the  Bloch  sphere  in  the  rotating  frame. 

One  can  use  these  single  qubit  rotations  not  only  as  part  of  the  collision,  but  also  for  initializing,  since  they  can 
bring  the  qubit  state  to  anywhere  on  the  Bloch  sphere.  As  already  discussed  in  section  ITvaI  the  ground  state  of 
a  coupled  PC  Qubit  system  is  not  in  general  the  product  of  single  qubit  ground  states.  Thus,  when  initializing  a 
qubit  via  Rz(0)  and  Rv{8 ),  one  is  not  starting  rotation  from  the  single  qubit  ground  state.  However,  since  the  ground 
state  is  very  close  to  a  product  of  single  qubit  ground  states,  this  difference  is  nearly  negligible.  In  figure  0  we  show 
the  effects  of  incorporating  this  error  into  the  algorithm  when  the  coupling  constant  is  taken  to  be  1/10  of  the  qubit 
resonant  frequencies,  a  rather  exaggerated  estimate  since  the  coupling  is  usually  much  smaller.  The  diffusion  constant 
is  decreased  by  this  approximation,  due  to  the  enhanced  population  in  the  1 11)  state  relative  to  the  |00)  state  from 
the  coupling. 

The  other  gate  operation  needed  to  form  a  universal  set  for  general  decomposition  is  coupled  free  evolution.  The 
first  thing  to  do  is  go  to  a  co-rotating  frame  where  one  can  have  a  coupled  Hamiltonian  only,  i.e. ,  no  single  qubit 
terms,  that  is  time- independent.  In  NMRQC  this  is  done  by  going  to  the  frame  where  both  qubits  are  rotated  around 
the  z  axis.  However,  since  our  coupling  does  not  commute  with  our  single  qubit  terms,  a  different  method  will  be 
used.  For  notational  convenience  only,  we  consider  the  case  where  both  qubits  are  biased  at  /  =  1/2,  where  our 
Hamiltonian  is 


H 


=  wxJl  +  will  +  J12L1/?1 


2t r 

h 


(ii) 


In  the  co-rotating  frame  where  both  qubits  are  rotating  around  the  x  axis,  one  has  the  Hamiltonian 


h  =  -Mil  i l 


■  i1!2] 

y  yi 


(12) 


as  long  as  wl  =  wi-  This  constraint  of  w\  =  wl  imposes  limitations  on  some  NMRQC  initialization  schemes  which 
use  frequency  selective  initialization. 

One  can  now  rewrite  the  unitary  collision  transformation  in  the  following  suggestive  way: 


Vswap  =  exp H^(<5/<j?  +  d^)]  exp 


(13) 
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The  first  term  is  just  free  evolution  in  the  co-rotating  frame.  The  second  term  can  be  written  as: 


expH-u^]  =  pi-r-alaWi-^RK--) 


(14) 


where  the  middle  term  can  be  written  as: 


exPh +  &ld'l)\Rl(7T) 


(15) 


Hence  one  can  perform  a  decomposition  of  the  collision  transformation  into  a  sequence  of  single  qubit  rotations  and 
coupled  free  evolution. 

In  summary,  we  have  shown  that  the  PC  Qubit  can  implement  the  unitary  transform  that  performs  collisions  in 
the  ID  FQLGA  for  the  diffusion  equation  by  a  single  and  coupled  qubit  evolution  decomposition.  The  single  qubit 
rotations  were  shown  to  be  feasible  for  qubit  initialization  as  well,  with  a  slight  approximation  due  to  the  coupled 
ground  state  that  is  not  a  product  state.  The  coupled  free  evolution  was  seen  to  require  identical  qubit  frequencies 
over  a  lattice  site,  making  initialization  a  bit  more  challenging. 


V.  CONCLUSIONS 

In  this  paper  we  have  shown  that  the  implementation  of  the  FQLGA  for  the  ID  diffusion  equation  is  feasible 
with  PC  Qubits.  We  began  by  considering  the  simplest  scheme  possible  using  the  PC  Qubit.  This  consisted  of  first 
initializing  the  qubits  while  keeping  them  in  their  ground  state,  and  then  performing  the  collision  by  quickly  changing 
their  flux  bias  points  and  then  performing  a  single  tt/2  pulse.  This  initialization  technique  could  prove  useful,  but 
the  way  we  have  implemented  the  collision  is  not  easily  generalized  to  other  collisions.  We  needed  to  develop  a  more 
general  collision  scheme,  and  then  see  how  we  could  initialize  in  conjunction  with  that  new  scheme. 

A  more  general  collision  transformation  was  then  discussed  by  decomposing  the  unitary  matrix  into  a  sequence  of 
single  qubit  rotations  and  coupled  free  evolution.  We  first  developed  single  qubit  rotations  for  the  PC  Qubit  that  could 
be  used  as  part  of  the  collision  decomposition  as  well  as  for  initializing  the  occupation  probabilities.  The  initialization 
was  considered  only  approximate  due  to  the  permanent  non-commuting  coupling  between  qubits.  For  the  coupled  free 
evolution  we  saw  that  transforming  to  a  rotating  frame  analogously  to  NMRQC  set  a  strong  but  feasible  constraint 
on  the  frequencies  of  our  qubits.  Ultimately  one  would  like  to  remove  the  constraint  of  equal  frequencies,  so  that 
frequency-selective  initialization  can  be  done  analogously  to  the  NMRQC  implementation,  alongside  the  very  general 
collision  scheme.  One  would  then  also  need  to  account  for  initialization  pulses  rotating  states  from  a  non-product 
ground  state. 
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